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Abstract 

We consider a generalization of the Jarzynski relation to the case where the 
system interacts with a bath for which the temperature is not kept constant 
but can vary during the transformation. We suggest to use this relation as 
a replacement to the thermodynamic perturbation method or the Bennett 
method for the estimation of the order-order surface tension by Monte Carlo 
simulations. To demonstrate the feasibility of the method, we present some 
numerical data for the 3D Ising model. 

1 Introduction 

In the last decade, important progresses have been made in the context of far- 
from-equilibrium statistical physics with the derivation of the so-called fluctuation 
theorems. Among these theorems, the Jarzynski relation has certainly become the 
most famous one and is now widely applied both in experimental and numerical 
studies not only in physics but in chemistry or biophysics as well. Remarkably, 
this relation relates thermodynamic equilibrium quantities to out-of-equilibrium 
averages over all possible histories of the system when submitted to a well-defined 
protocol. The system is initially prepared at thermal equilibrium and then driven 
out-of-equilibrium by varying a control parameter h from say hi to /12. The work 
W is recorded for each experiment. The Jarzynski relation states that 

e -0AF = {e -0 W) (1) 

where AF is the free-energy difference AF — F{1%2) — F(hi) between the two equi- 
librium states at values hi and hi of the control parameter. In his first paper on the 
subject, Jarzynski gave a derivation of this relation in the case of a set of particles 
evolving according to the laws of Newtonian mechanics [1] . The derivation depends 
crucially on the assumption that the system is initially thermalized with a heat bath 
at temperature T but then isolated from it while the work is exerted. No exchange 
of heat with the environment is possible. A change of the temperature of the bath 
has thus no consequence. In a second paper [2] , Jarzynski showed that the relation 
(1) may be established in the case of a Markovian dynamics too. In this case, the 
interaction with the bath is encoded in the transition rates. The demonstration 
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does not require anymore the insulation of the system from the bath while the work 
is exerted. Exchange of heat with the environment is properly taken into account 
during the whole process. 

It is a simple exercise to generalize the Jarzynski relation to the case where the 
temperature of the heat bath changes with time. In the special case where no work 
is exerted on the system, one obtains (24): 

e -A(f3F) = , e - J E(t)f)du ^ 

where E(t) is the total energy of the system at time t. The derivation can be found 
in the appendix. As noted by Crooks in his dissertation [3] "it is possible to view 
this change as an additional perturbation that exerts a sort of entropic work on the 
system" . Like any other Jarzynski relation, Eq. (2) is also a special case of the 
annealed importance sampling method [4] . As far as we know, the case of a varying 
temperature has not attracted any attention. The aim of this paper is to show that 
a Jarzynski relation generalized to a temperature change may also have some inter- 
esting applications and deserves interest. As an example, we propose a new way to 
estimate by Monte Carlo simulations the order-order surface tension. This method 
is tested for the three-dimensional Ising model. The surface tension is related to the 
remaining free energy of a system containing two phases in coexistence separated 
by an interface when the contribution of these bulk phases has been removed. The 
Jarzynski relation offers the possibility to calculate these free energies in a Monte 
Carlo simulation. 

The precise numerical determination of the surface tension has attracted a con- 
tinuous interest for several decades. It is indeed an important quantity in nucleation 
theory and may be used to distinguish between a first-order phase transition and a 
continuous one. In the latter case, the surface tension is expected to vanish at the 
critical point. In the context of lattice spin models, the interface tension permitted 
to show for example the existence of a randomness-induced continuous transition 
in the 2D eight-state Potts model [5] or that the phase transition of the 3D 3-state 
Potts model is weakly first-order [6]. In chemistry, one is interested in the surface 
tension in fluids or colloids and the 3D Ising model is commonly used as a toy model 
for the description of the liquid- vapor phase transition [7, 8] . In the context of high- 
energy physics, the confincmcnt-dcconfinemcnt phase transition between hadronic 
matter and quark-gluon plasma is a very active research field. Numerous Monte 
Carlo simulations have been devoted to the estimation of the interface tension of 
the pure SU(3) gauge model [9, 10, 11, 12, 13, 14]. The Ising model is also of in- 
terest in this context since duality maps the 3D Zi Yang-Mills gauge theory onto 
the 3D Ising model [15]. Wilson loop and Polyakov loop correlators are then re- 
lated to interfaces in the Ising model. In all theses situations where the surface 
tension is of interest, numerical accuracy has been constantly improved not only by 
the exponential growth of computer speed but mainly by progresses in algorithms 
(multicanonical simulations [16], flat histogram methods [8],. . .) and in the protocol 
set up to have access to the surface tension (thermodynamic integration, Binder's 
histogram method [17], snake algorithm [18, 19, 20, 21], . . .). 

This paper is organized as follows: in the first section, we expose the method 
used to estimate numerically the order-order surface tension. The derivation of 
the Jarzynski relation is given in the appendix. Numerical results for the 3D Ising 
model are presented in the second section. Conclusions and discussion follow. 
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2 Numerical estimation of the surface tension 



Consider the classical ferromagnetic Ising model denned by the usual Hamilto- 
nian [22] 

Wising = — J ^ CjCTj, Ui = ±1 (3) 

(i,3) 

where the sum extends over nearest neighbors only. In order to favor the appearance 
of an interface in the system, we first impose anti-periodic boundary conditions in 
one direction. The free energy F(T) of the system at a temperature T < T c may 
be decomposed as 

F(T) = F Ferro (T) + F S (T) - k B T\nL (4) 

where F-p „ is the contribution of the two ferromagnetic phases on both sides of 
the interface, F s the free energy of the interface and k B lnL the entropy associated 
to the degeneracy of the position of the interface. One can restrict the analysis to 
the case where only one interface appears in the system. Indeed, the probability 
for n interfaces behaving as e ~ /3 " cri , the contribution of spin configurations with 
more than one interface can be neglected for sufficiently large lattice sizes L. The 
order-order surface tension is defined as 

a(T)= lim (5) 

To first calculate the free energy F(T), we will estimate numerically the differ- 
ence A([3F) = -^jr — kB jr . The temperature Tj <C T c is chosen such that the 
thermal contribution to the free energy F{T{) can be neglected, i.e. F(Ti) — 
E Q - fc B Tln2 - k B T\nL where E Q = -J(dL d - 2L d ~ 1 ) is the energy of the spin 
configuration with a flat interface, ks In 2 the entropy associated to the degeneracy 
of the ground state corresponding to Ui — > — at and k B InL the entropy associated 
to the position of the interface. 

The usual method to compute A(/3F) is based on the decomposition 

e -A(/3F) = = Z Pn + i Z f j N Zp 2 ^ 

where {/3\, . . . , f3jy + i} is a set of inverse temperatures interpolating between /3j = (5\ 
and (3 = (3n+i- When (3k and fik+i are sufficiently close, the ratio Zp k+1 /Zp k can 
be estimated during a single Monte Carlo simulation at the inverse temperature (3k 
by thermodynamic perturbation: 

ip±L = {e -(h + i-h)E )h (?) 
Z Pk 

or by more elaborated procedures like the Bennett method [23] . The estimation of 
A((3F) requires N Monte Carlo simulations at the inverse temperatures (3\, . . . (3n 
with for each of them a sufficiently large number n cxp . of Monte Carlo steps in order 
to get accurate averages. 

The Jarzynski relation (23) offers the possibility to estimate the quantities 
Z/3 k+1 /Zfj k . The system is initially thermalized at the temperature Tk = l/k B Pk- 
First the initial spin configuration is stored. Then niter. MCS are performed with 

a temperature increasing from Tk to Tk+\. The quantity e "Her. Z^*=o 

is calculated and accumulated. The initial spin configuration is restored and a 

few additional MCS at temperature Tk are performed to generate a new initial 
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spin configuration uncorrelated with the previous one. The whole procedure is 

repeated n cxp . times. Finally the average (e "iter. Z^ t =o 1 gives e~ " ' 
where A(/3F) = j3k+iF([3k+i) — (3kF(j3k) according to the Jarzynski relation (23). 

Both methods lead to a numerical estimate of F(T). The last step is now to 
calculate F FeTTO (T), i.e. the contribution of the two ferromagnetic phases. To that 
purpose, we make a second Monte Carlo simulation with periodic boundary con- 
ditions. The method presented above is applied again to estimate the difference 
A(/3F Forro ) = Ff — j T) - %^p- The approximation F Fcrro (T,) = E - k B T\n2 
where E — —JdL d gives finally Ff ciio (T). We can now use (5) to estimate the 
surface tension. 

The algorithm based on the Jarzynski relation depends on three parameters: N 
the number of temperature intervals in which [/3^ ; /3] is divided, n^er. the number 
of MCS bringing the temperature from Tk to Tk+\ and n oxp . the number of mea- 
surements of e _A ^S E performed. The thermodynamic perturbation method (7) 
corresponds to the special = 1. The use of the Jarzynski relation may 

improve the convergence and save computer time by an optimal choice of the set 
(N, niter., n exp .) of parameters. Increasing niter, allows to decrease N and to some 
extend n oxp .. However, we have no recipe to determine these optimal parameters 
and moreover, they probably depend a lot on the model and on the Monte Carlo 
dynamics chosen. 

3 Application to the 3D Ising model 

The method presented above is applied to the 3D Ising model. We considered 
the "extreme" case N = 1, i.e. the temperature is not increased step by step in 
different Monte Carlo simulations but in a single simulation from Tj to T. This 
choice is certainly not the optimal one but it demonstrates that the Jarzynski re- 
lation can be applied even when the temperature is changed by a very large amount. 

In one of the three directions, the boundary conditions are first anti-periodic 
and then periodic. In the two other directions, we are free to choose any bound- 
ary conditions. We considered both periodic (PBC) and free boundary conditions 
(FBC). This choice is motivated by the fact that only periodic boundary conditions 
are usually considered in the literature but in the 2D Ising model, one can easily 
check that finite-size corrections to the surface tension are smaller with FBC. This 
can be understood by considering the equivalent one-dimensional Solid-on-Solid. 
With FBC, the steps made by the interface are independent and the free-energy 
may be written as the sum: 

F s = -k B TL\n[ ^ e-W+W)" 

n— — oc 

= 2JZ. + fc B T£lntanh/3J (8) 

which leads to Onsager's exact result a = 2 J + fcsTlntanh/3J for any lattice width 
L. When the lattice size perpendicular to the interface is finite, corrections arise 
because the sum of (8) becomes bounded. PBC impose a non-local condition that 
results into additional finite-size corrections. 

We considered cubic lattices with sizes L = 4, 6, 8, 12, 16, 24, 32,48. The initial 
inverse temperature is = 2, to be compared with the critical inverse tempera- 
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ture f3 c ~ 0.22165 . . . [24]. The spins are initialized in the state <7j = +1 if they 
lie in the upper half of the system and in the state <7j = — 1 otherwise, i.e. we 
start with a flat interface separating two ferromagnetic phases at saturation. We 
made 500 x L MCS to equilibrate the system at this temperature. This choice 
is a very safe bet since in practise very few spin flips occur during these itera- 
tions 1 . The inverse temperature (3 is then decreased up to the final value in the 
range [0.24; 0.80]. The final temperature is reached after niter. = 250 x L MCS. 
The experiment is repeated n cxp . = 1000 times. Before each experiment, the spin 
configuration is stored. After the experiment, it is restored and 25 x L additional 
MCS are performed to generate a new non-correlated spin configuration for the next 

experiment. The statistical error on the average is estimated as the 

standard deviation as expected from the central limit theorem. However systematic 
deviations may occur due to the fact that the average is dominated by rare events. 
This effect is important when the variation of the temperature is fast, i.e. when 
the transformation is strongly irreversible. Moreover, the bias induced is difficult 
to quantify. As shown by Gore et al. [25, 26], "the number of realizations for con- 
vergence grows exponentially in the average dissipated work" . To overcome this 
difficulty, we took advantage of the fact that the detailed balance condition (16) is 
the only assumption made on the transition rates in the derivation of the Jarzynski 
relation. Any Monte Carlo algorithm satisfying this condition can thus be used. 
To make the transformation as reversible as possible, we used cluster algorithms: 
the Hasenbush-Meyer algorithm [27] for anti-periodic boundary conditions and the 
Swendsen-Wang algorithm [28] for periodic boundary conditions. 
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Figure 1: Surface tension a versus the number of iterations niter, for three different 
temperatures: f3 = 0.24, (3 = 0.40 and (3 — 0.56 from bottom to top. The lattice 
size is L = 16. Error bars on the data correspond to statistical errors. 



The surface tension is calculated according to the method presented in section 2. 
Estimates of the surface tension for an increasing number of iterations niter, are plot- 
ted on figure 1. The number of experiments is kept fixed to the value n CX p. = 1000 

lr The probability that a spin be flipped during a Monte Carlo step is e~ 4d ''*' / where d is the 
dimension of the lattice. At the inverse temperature ft J = 2 and dimension d = 3, this probability 
is as small as 10~ 14 . 
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which means that rare events with a probability smaller than l/n cxp . = 10~ 3 are 
not correctly sampled. The figure shows two regimes. When the number of itera- 
tions is large n n it or > n°^ T , the estimates display a plateau. The amplitude of the 
fluctuations around the value of this plateau is of the same order of magnitude than 
the error bars. In contradistinction, the estimates display a systematic deviation 
when the number of iterations is small n n i te r < ^° t er • The data indicate a deviation 
proportional to In - ^ . The contribution of the rare events are larger than the 
statistical errors and cannot be neglected. The cross-over between the two regimes 
depends on the temperature. Surprisingly, it does not depend much on the lattice 
size. 

We have also tried to estimate the surface tension using the width a 2 Vs = (Wj) — 
(Ws) 2 of the probability distribution of the entropic work. These distributions are 
indeed close to Gaussian distributions so the Jarzynski relation reads 

{e -w s) = _J_ [ e-^-^-^'^sdWs^ A((3F) = (W S )-^- (9) 

The estimates of the surface tension are compatible with those obtained using the 
Jarzynski relation (2) but the error bars are larger. As a consequence, we will not 
consider this estimate in the following. 

Our numerical data of the surface tension are presented on figures 2 (PBC) and 3 
(FBC). The few points with large error bars on figure 3 correspond to small systems 
at temperatures higher than the finite-size "critical" temperature T C (L). The system 
is thus already in the paramagnetic phase and the interface disappears. No order- 
order surface tension can then be given. In contradistinction to the two-dimensional 
case, the surface tension displays larger finite-size corrections with FBC. The data 
have been extrapolated to the thermodynamic limit using the ansatz: 

a(L) = a(oo) + ^ + § (10) 

for PBC and 

(7 (L)=a(oo) + | + ^ (11) 

for FBC. The error bars on the value ct(oo) take into account both the errors on 
the data and due to the fit. The figure 4 gives an example of this extrapolation 
for three different temperatures. While the whole range of lattice sizes L are well 
reproduced by the ansatz (10) for PBC, the two smallest lattice sizes (L = 4 and 
L = 6) cannot be taken into account by (11) for FBC and have been discarded. 
The extrapolated value cr(oo) of the surface tension can be compared in figures (2) 
and (3) to the 17th-order low-temperature expansion (see [29] for the coefficients) 
and to large-scale Monte Carlo simulations [29] . Our data are perfectly compatible 
to the latter apart from a small disagreement for FBC at temperatures close to the 
critical point. The too simple ansatz (11) may be responsible for this disagreement. 
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Figure 2: Numerical estimates of the surface tension of the 3D Ising model with 
respect to the inverse temperature for periodic boundary conditions. 
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Figure 3: Numerical estimates of the surface tension of the 3D Ising model with 
respect to the inverse temperature for free boundary conditions. 



4 Discussion and conclusions 

We have discussed a Jarzynski relation generalized to temperature changes. Its 
derivation is limited to Markovian dynamics since the interaction with the bath is 
properly taken into account only in this case. We have then presented a modified 
version of the usual algorithm allowing for the computation of the surface tension by 
Monte Carlo simulations where the Jarzynski relation replaces thermodynamic per- 
turbation. The approach has been tested in the case of the 3D Ising model. Despite 
the common belief that rare events prevent the application of Jarzynski relation, 
the method is efficient because the transformation does not bring the system far 
from equilibrium. This is possible because the method uses initial and final states 
not very different. Basically, the final state display thermal fluctuations while the 
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Figure 4: Extrapolation of the surface tension for three different temperatures. The 
solid curves correspond to periodic boundary conditions and the dashed ones to free 
boundary conditions. 



initial state is almost the ground state. But in both cases, two ferromagnetic phases 
are separated by an interface. In contradistinction to other algorithms, the inter- 
face is not created during the transformation. Moreover, we use cluster algorithms. 
Their relaxation times are much smaller than the usual Metropolis algorithm and 
allows the final state to be close to equilibrium. The systematic deviation due to 
an insufficient sampling of the rare events that dominates the average {e~ Ws ) has 
been shown to be negligible compared to statistical errors for a reasonable number 
of iterations. Our aim was only to prove the usefulness of the Jarzynski relation 
despite the fact that rare events play an important role and not to add a 6th or 7th 
new digit to the estimates of the surface tension that can be found in the literature. 
The computational effort we devoted to the calculation is much smaller (around 
100 times) than that of rcf. [29] for instance. For (3 — 0.45 and L = 32, our esti- 
mate of the surface tension is 1.8316(4) while Hasenbusch et al. obtained a much 
accurate estimate: 1.831368(5). However, they only calculated numerically the dif- 
ference F{(3 = 0.45) — F((3i = 0.515) and then used the 17th-order low-temperature 
expansion estimate of F((3 = 0.515) to obtain the surface tension. We calculated 
numerically the difference F(/3 — 0.45) — F(fli = 2) and estimated F(/3i — 2) as 
the T = ground-state free energy. By starting with a much temperature, we do 
not need the 17th-order low-temperature expansion to calculate the surface tension. 

The main drawback of our approach is that the convergence of estimates of the 
free energy obtained by the Jarzynski relation is not well understood yet. We can- 
not give any recipe to adjust the three parameters N, ni ter ., «cx P . to obtain the best 
convergence. Research efforts in this direction are highly desirable. 

One can imagine other protocols based on the Jarzynski relation to estimate the 
surface tension. One can start first with a strong homogeneous magnetic field on 
the up and down boundaries to force the system in the ferromagnetic state. Then 
one of these magnetic fields is reversed to favor the appearance of the interface. By 
recording the work while the magnetic field is reversed, the free energy could be 
estimated using the more usual Jarzynski relation (1). However, the convergence 
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may not be as good because an interface has to be created during the magnetic field 
reversal. In the approach based on variation of the temperature, this interface is 
already present in the initial state. 

The relation (22) may be useful even in the case of the measurement of the 
difference of the free energies of two equilibrium states at the same temperature. 
If these two states are separated by a free energy barrier, the relation (22) allows 
first to increase the temperature so that the system can pass the barrier and then 
to decrease the temperature down to its initial value. 
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Appendix: derivation of the Jarzynski relation 

We shall derive the Jarzynski relation along the lines of the derivation given in the 
appendix of [2]. Let us denote a the microscopic states of the system, i.e. the spin 
configurations in the case of the Ising model. We shall assume that the dynamics 
of the system is Markovian and is governed by the master equation 

p(a, t + l) = J2 P(^> t)W mMt) {5' -> a) (12) 

where p(a,t) is the probability to find the system in the state a at time t and 
Wf3( t ).h(t){v —* is the transition rate per time step from the state a to the state 
a'. The transition rates take into account the interaction of the system with an 
heat bath in an effective manner. They depend on the temperature T = l/fcg/3 
of the bath and on an external parameter h, for instance the magnetic field. The 
master equation is equivalent to the Bayes relation which means that the transition 
rates are conditional probabilities: Wpn^M^a — ► <?') = p(a',t+ l\a,t). As a 
consequence, the transition rates satisfy the condition 

J2Wp, h (a ^ a') = I (13) 

When (3 and h are kept constant, the system is expected to evolve towards the 
Boltzmann equilibrium distribution 

^ h {S) = J_ e -/W)-fcM (<? )) (14) 

Zp. h 

where M is the extensive quantity associated to the control parameter h (for in- 
stance the magnetization when h is the magnetic field) . The equilibrium distribution 
satisfies the stationarity condition 

*p,h{3) = ]T n , h (a')W 0!h (a' -> 3) (15) 

3' 

For convenience one imposes the more restrictive condition of detailed balance 

*f},h(Z)W fl , h {ff -> <?') - n^hia^W^hia' -> a) (16) 
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that leads to (15) when using the condition (13). 

At time t = ti the system is thermalized with a bath at the temperature 
l/ksPiti) and external parameter h(ti), i.e. p(a,U) = T^fjUAMu)^)- The tem- 
perature and the external parameter are then varied and the quantity 

(x( ti , tf )) = (e -E;- l (A t ^(^))-A t(W M(,( t )) )) (17) 

is measured. The brackets (. . .) denotes an average over all possible histories of the 
system between ti and tf. We have introduced the notations A t f3 — (3(t + 1) — (3(t) 
and A t (/3h) = 0(t + l)h{t + 1) - (3(t)h(t). Following Jarzynski [2], each time 
step consists into two substeps. First the temperature and the external param- 
eter are changed by an amount A t (3 and A t h. The latter induces a work W — 
—A t hM. In the second substep, the master equation is iterated once so that the 
system makes a transition to a new state a(t + 1) with the conditional probability 
Wp( t +i),h(t+i){o{t) -» a(t + 1)). The energy change E(a(t + 1)) - E(a(t)) corre- 
sponds to the heat exchanged with the heat bath. Rewriting X(ti,tf) using the 
equilibrium probability (14) as 

t f l 

X ( t . t j = TT Z fHt+i)Mt+i) 1T (i{t+i).h{t+^^ t )) ^ 
t =7. 2/3(t),h(t)7r0(t),/i(t)(0W) 

the average (17) reads 

1 1 — l 

{X{U,t f )) = 2^ n f3(u)Mu)( a W) 11 5 T^TP\\ 

xW 0{t+1)Mt+1) (a(t) ^a(t + I))}l9) 



All but two partition functions cancel. Using the detailed balance condition (16), 
one obtains 



{X{ti,tf)) = — 2^ ^(tOMnn^W 11 ~ { rritW ^ 



x W [j{t+1)Mt+1) {a{t + 1) -> 
Now all but one equilibrium probabilities cancel and one gets 

- / (t/) ' Mt/) £ *m)Mt f Mtf)) II ^(t+i),Mt+i)(^(*+ 1 ) ->*(*)) 

^/9(t*),h(t«) {s{t)} t=u 

(21) 

Finally the sum over histories of the system is easily shown to be equal to 1 using 
equation (15). Introducing the definition of the free-energy F(f3, h) = — /3 _1 \nZp^, 
one obtains the Jarzynski relation 

^-£^7/(A«/»E(*(t))-At(/"0^(*(t)))j = e -[/s(t f )F(0(t,)Mtf))-0(ti)F(f3(u) M ti))] ( 22 ) 
In this paper, we restrict ourselves to the case h = 0, i.e. 

( e - EtC A t ^B(ff(t))j = e _[ /3 ( t/ )F(/3(t / ))-/3(t l )F(/3(t I ))] ^3) 
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The relation (22) is appropriate for Monte Carlo simulations since the time is 
discrete. The demonstration given by Jarzynski in the case of a Markovian process 
in continuous time can also be extended to a varying temperature [2] . It leads to 

(e - m [ms(t))-PhM { s md t ) = e - [mF{mMm if (24) 

where U(a(t)) = E{a(t)) — h{t)M{a{t)) is the total energy of the system in the 
state a(t). 
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